Dynamic Localization of Interacting Particles in an Anharmonic Potential 



M. Herrera/'0T. M. Antonsen,!'^ E. Ott,^'^ and S. Fishman'"^ 

^Dept. of Physics and IREAP, University of Maryland, College Park, Maryland 20742, USA 
^Dept. of Electrical and Computer Engineering, University of Maryland, College Park, Maryland 20742, USA 
Physics Dept., Technion- Israel Institute of Technology, Haifa 32000, Israel 
(Dated: December 13, 2012) 

We investigate the effect of anharmonicity and interactions on the dynamics of an initially Gaus- 
sian wavepacket in a weakly anharmonic potential. We note that depending on the strength and 
sign of interactions and anharmonicity, the quantum state can be either localized or delocalized 
in the potential. We formulate a classical model of this phenomenon and compare it to quantum 
simulations done for a self consistent potential given by the Gross-Pitaevskii Equation. 
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Atoms inside anharmonic traps display a variety of 
phenomena including wavepacket spreading due to de- 
phasing, and quantum revivals pQ [5] . Additionally, sys- 
tem behavior is also strongly influenced by particle- 
particle interactions: for example, theoretical [3j and ex- 
perimental [4] studies have demonstrated that interac- 
tions alter the frequency of collective modes of a Bose- 
Einstein Condensate (BEC). Recently, systems with both 
anharmonicity and interactions have been studied. The 
effects of interactions and anharmonicity on the stability 
of stationary states O |6], collective motion [TlfTO] and 
dynamics of coherent states [TT] of BECs have been in- 
vestigated analytically and numerically. Additionally, ex- 
perimental work has also studied the dynamics of BECs 
in the presence of anharmonicity [T^HIl] . 

Here we consider the quantum evolution of an initially 
Gaussian wavepacket in a one-dimensional trap with both 
interactions and anharmonicity, using the description 



y 



(1) 



where Va = Px'^/'i, and Vint = m|V'P- This can be consid- 
ered as the mean field approximation for the condensate 
wavefunction 'ijj{x,t) of a BEC with interactions, and is 
known as the Gross-Pitaevskii (CP) Equation [T5l417j . 
Here /3 quantifies the anharmonicity of the trap, and u is 
a nonlinear interaction strength, which quantifies repul- 
sive (u > 0) or attractive {u < 0) interactions. We con- 
sider the regime where Va, Vint are smaller than the har- 
monic oscillator Hamiltonian, and adopt harmonic units, 
X ~ X I (^^/hjrrvj^o) , t = ujQt, where x and t are the un- 
normalized units, wq is the frequency of the harmonic 
oscillator, m is the mass, and the position and momen- 
tum operators, x and p satisfy [x,p] = i, with the nor- 
malization J dx\'ip\'^ = 1. We numerically solve Eq. ([l]) 
using a split-step Fourier method [TS]. One possible ex- 
perimental realization is the use of a Feshbach resonance 
[19H22] in order to sample small repulsive and attractive 
interaction strengths. Consider an initial, approximately 
Gaussian wave-packet centered at x = Xc, with symmet- 
ric widths in x and p, — Ap = •y/l/2. This state 



may be formed by taking the confined ground state at 
the center of the potential and discontinuously shifting 
the potential to the left by an amount Xc- In the case 
where the potential is formed by the interference of two 
laser beams, this shift might be accomplished by chang- 
ing the relative phase of the two beams, for example as 

in Eg. 

In the absence of interactions, a wavepacket that over- 
laps many energy eigenstates of an anharmonic trap will 
dephase due to the nonuniform spacing of the relevant 
energy eigenstates. A classical interpretation is that in 
the phase space, an initial Gaussian distribution subtends 
various action tori, each with a different frequency. Thus, 
different parts of the distribution function rotate in the 
phase space at different rates, leading to spreading of the 
initial distribution function. 

In the presence of interactions, however, we have nu- 
merically observed what we term a "dynamical localiza- 
tion" phenomenon in the time evolution of the position 
density of the wavepacket. For some choices of inter- 
action strength u and anharmonicity /3, the wavepacket 
does not spread over the trap, while for others it does. 
This is illustrated in Fig. [T where panels (a) to (c) il- 
lustrate the position density -01^7 obtained from solution 
of Eq. 0, at late times {t = 1999) for /3 = 1.89 x 10"'' 
and three values of interaction strength, one attractive 
(a) and two repulsive (b and c). Only in the weakly re- 
pulsive case (b) does the wave function spread through- 
out the trap as it would in the absence of interactions. 
Figure displays the value of cTxCTp, where ax {op) 
is the standard deviation of position (momentum), av- 
eraged from t = 1899 to t — 1999, for various choices 
in u and /3, and Xc — —8. A localized state will have 
smaller values c^CTp, shown as dark regions while a de- 
localized state will have larger values, shown as bright 
regions. The solid (blue) straight lines in Fig. [ijd) are 
results of an approximate theory (given below) for the 
parameter boundaries separating localized and delocal- 
ized behavior. It is important to note that, though one 
might intuitively expect attractive interactions {u < 0) 
to always lead to localization. Fig. [l] demonstrates that 
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this is not so. In fact, for a given value of u, both local- 
ized and delocalized behavior is possible, and reversing 
the sign of both u and /3 leads to the same behavior. To 




-10 10 -10 10 




-10 10 -0.5 0.5 



X u 

FIG. 1. (Color Online) (a)-(c): The position density lipl^ 
at t = 1999 for /3 = 1.89 x 10"" with different values of u, 
(a) u = -0.2586, (b) u = 0.1552, and (c) u = 0.50. Note 
that increasing repulsive interaction strength can lead to lo- 
calization, (d): Values of a^Up for different values of u and 
j3. Bright regions indicate delocalized behavior (large a^o-p), 
dark regions indicate localized behavior (small cJxOp). The 
solid lines (blue) are theoretical predictions separating regions 
of poor initial confinement and strong initial confinement. 

gain insight into the results represented in Fig. [TJ we in- 
troduce a classical model for the dynamical localization 
effect, which is motivated by our initial choice \xc\ — 8, 
in terms of a classical Hamiltonian. We note for the case 
of Fig. [T] where the state before the shift of the poten- 
tial is approximated by the ground state of the harmonic 
oscillator, that, after the shift of the potential, one can 
show that the state overlaps on the order of \xc\ energy 
levels. Hence, if \xc\ substantially exceeds one, we ex- 
pect that a classical treatment will be relevant. Thus we 
start by considering the classical phase space evolution of 
the distribution function f{x,v,t) of a one-dimensional 
harmonic oscillator with Hamiltonian Hq = + a;^/2. 
Suppose the distribution function is initially Gaussian (in 
analogy to the initial Gaussian wavepacket in the quan- 
tum case) and is displaced away from the bottom of the 
potential so that it is centered about a point in phase 
space (xc,0). In the absence of interactions and anhar- 
monicity, this distribution coherently orbits in the phase 
space, since the frequency is independent of the action. 
We introduce as small perturbations an anharmonic po- 
tential Va = /?a;^/4, as well as a self interaction potential 
Vint = un{x,t) where n{x,t) = J dvf{x,v,t) is the den- 



sity. The Hamiltonian is then H = Hq + eVa + (Vint, 
where e is an order counting parameter. Since Va and 
Vint are small perturbations to the harmonic oscillator 
Hamiltonian Hq, the center of the phase space distribu- 
tion continues to oscillate in the potential at a frequency 

close to the unperturbed value of one. Changing to 
the action-angle variables of the unperturbed oscillator, 
J and 9, and moving to a rotating frame (p — 9 — fit, 

the equations of motion become 

dJ/dt = -dH/d(t>, d(f>/dt ^ dH/dJ 
H = J{1 - fl) + epj^ sin^(</> + VLt) 

+ ey„t(\/2Jsin(0 + l]i),i), (2) 

where p = V^J cos{9) and x = \/2Jsm{9). Recalling 
that Va and Vint are small perturbations to the harmonic 
oscillator, we note that there are two different time scales: 
a fast time corresponding to the center of mass motion of 
the distribution function, and a slow time over which the 
shape of the distribution function evolves. We implement 
a separation of time scales t to, and ti — eto, J — 
Jq + cJi and </> = ^o + e^^i- We also expand the oscillation 
frequency of the center of the distribution function as 
~ l+e6il. Averaging the resultant equations of motion 
over the period of the fast time scale, and requiring that 
{dJi/dtQ) = {dcpi/dtQ) = (i.e. the correction terms do 
not grow secularly as functions of to) yields equations for 
the slow time evolution: 

dJo/dti = -d {H)/d(t)Q, d<j>Q/dti = d {H)/dJo 

(if) = ^/3Jo Qjo- Jc) (3) 

where Jc = The evolution of the slow-time trajec- 

tories (or the trajectories averaged over the short period) 
approximately follow contours of constant {H). We have 
determined the value of 5^ by requiring that = 
since we are in a moving frame such that the center 
of the distribution function is stationary, and assuming 
no contribution to the frequency from (Vi„t). In order 
to average the interaction potential, we are required to 
average the density over an oscillation period. Recall 
that the initial phase space distribution is Gaussian in x 
and p, centered about some x = Xc, p = Q with width 
Aj. = Ap = A = a/1/2. We now evaluate {Vint) for the 
early time evolution of the phase space distribution, i.e., 
when that distribution is still approximately Gaussian. 
Noting that the Gaussian distribution does not change 
appreciably over one oscillation period, allows us to av- 
erage Vint, leading to 

(Vi„t) = - — \ dioexp[-(\/2Jsin((/)o + to) 

-y2Xsin(to))']. (4) 

Numerically calculating the average of the potential, 
we can then estimate (H) . Figuresl2ja)-(d) plot contours 
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FIG. 2. (Color Online) Contours of {H) (black curves) for 
various values of interaction strength u (a)=-0.20, (b)=-0.02, 
(c)=0.10, (d)=0.5 with /3 = 2 x 10"*. Separatrices are shown 
in bold, separating free streaming and confined trajectories. 
In red are numerically calculated trajectories under the in- 
fluence of an external potential of the form in Eq. Q. The 
trajectories consists of short time (fast oscillations) and long 
time behavior. 



of (H) in black for different choices of interaction strength 
u and anharmonicity /?, with = 32, \xc\ = 8. 

There is a difference in the nature of the phase space 
when u and /3 have opposite signs, Figs. W[a.) and (b), 
and when u and /3 have the same sign, Figs.l2lc) and (d). 
From Eqs. ([3| and Q, we see that reversing the sign of 
both u and /3 leaves phase space unchanged. For u and 
P of opposite signs, a fixed point is present at ( J, </>) — 
(Jc,0), and a separatrix (shown in bold) separates free 
streaming trajectories from trajectories orbiting the fixed 
point. In contrast, for u and f3 of the same sign, while an 
elliptic point is still present at (J, (f)) = { Jc, 0), additional 
hyperbolic points typically appear on the = axis. A 
separatrix (bold) again separates trajectories which orbit 
about the fixed point from those that do not, and, the 
trajectories immediately outside the separatrix are swept 
away from the region of the fixed point to values of higher 
J. Note that the size of the region of confined trajectories 
increases for larger values of |m//3|. Plotted in color in 
Fig. [2] are numerically calculated trajectories with Vint 
given by Eq. [2j As can be seen, the trajectories have 
both a short time motion (fast oscillations) and long time 
motion in the phase space, with the long time motion 
given approximately by the contours of {H) . 

We interpret the formation of a classical elliptic re- 
gion around (J, <^) = (Jc,0) as arising from a nonlinear 
resonance P5! between the integrable motion of the an- 



harmonic well and a periodic driving term. This has been 
studied quantum mechanically, for example, in the rovi- 
brational modes of diatomic molecules j25) and in the 
dynamics of Rydberg atoms [57] . One important dif- 
ference, is that in this example, the driving is not due 
to an external field. Rather, the condensate/ wavepacket 
drives itself on resonance because Vint = un{x,t). As 
we discuss below, this feedback between the density and 
Vint can lead to dynamics which can either be localizing 
or delocalizing. 
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FIG. 3. (Color Online) The Husimi functions at t = (a)- 
(b) , and at t — 1999 (c)-(d) under the influence of the CP 
equation. In (a) and (b), the classical separatrix from the 
early time model is plotted in bold. In (a) u = 0.5, /3 = 
2 X 10"**, the initial condition is well contained inside the 
classical separatrix found with the early time model . In (b) 
u — O.l, 13 = 2 X 10""*, however, the initial condition extends 
well beyond the classical separatrix. (c) demonstrates that 
the initially well contained wavepacket continues to remain 
localized at late times, while (d) demonstrates that the poorly 
contained wavepacket has spread in the potential (color scale 
narrowed to show detail). (In (c) we have applied a small 
shift Aff) = 0.31 to recenter the Husimi function on <jf> = 0. ) 



Figure |2] is constructed using a prescribed oscillating 
Gaussian potential. We now replace it with the self con- 
sistent interaction potential Vint — un{x,t). The con- 
tours of (H) found from Eqs. ([3| and (|4| will only 
describe classical trajectories for early times unless the 
wavefunction retains its localized character. To this end 
we consider the size of the confining region compared to 
the extent of the Wigner function (the quantum analog 
to the phase space distribution, which is equal to the 
Gaussian distribution function at t = 0). 

Building upon insight gained from the classical model, 
it is reasonable to expect that if an initial Wigner func- 
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tion is well confined, that is, enclosed in the classical sep- 
aratrix denoting orbiting trajectories, the Wigner func- 
tion will remain localized in the phase space at later 
times. Similarly, if the Wigner function, is not well con- 
fined in the separatrix, then at later times it will be- 
come delocalized in the phase space. Qualitatively, as 
the wavepacket leaves the confining region, the strength 
of Vint decreases. This altering of the potential leads 
to time dependent behavior with more of the wavepacket 
leaving the confining region, and coexistence of a delocal- 
ized component of the state, and a localized component 
(most easily seen in the Husimi function representation) 
at lower values of action. 

For example, in Figs, [sjja) and [sjjb) we plot the ini- 
tial Husimi function Q( J, (j>) (equivalent to the Wigner 
function smoothed over a Gaussian function in order 
to assure non-negative values for all times), in the (J, 
(j)) coordinates at < = for various values of u with 
/3 — 2 X 10~*. Plotted in bold is the classical separa- 
trix from the Gaussian model described previously, for 
the figures' respective values of u and /3; in Fig. [Sj^a) 
u = 0.5 and ;3 = 2 X 10-^ while in Fig. ^h) u = O.l 
and /3 = 2 X 10~^. We note that in Fig. [sjja), the Husimi 
function is still relatively localized inside the classical sep- 
aratrix at late times. Fig. [sj^c). However, in Fig. [sjb), 
where the initial Husimi extends beyond the classical sep- 
aratrix, the Husimi function it not as well localized in- 
side the classical separatrix at late times Fig. [sf^d). The 
quantity F appearing in Figs, [sjja) andjsj^b) is the frac- 
tion of the initial Wigner function which lies inside the 
confining separatrix, and reflects the difference in initial 
confinement between the two cases. 

We test the separatrix expectation numerically by 
choosing different values of u and /?, and calculating the 
quantum evolution of an initially Gaussian wavepacket 
centered at {x,p) = (—8,0). Recall that the classical 
model demonstrates that the size of the confining re- 
gion is dependent on the sizes and signs u and (3. Using 
our classical model, we find for given values of u and /3 
the fraction of the initial Wigner function that is con- 
tained in the confining separatrix, F{u,f3). Figure [l|d) 
plots in color (bold lines) contours of F{u, f3) correspond- 
ing to = 0.9. Note that regions yielding poor initial, 
classical confinement, F < 0.9 correspond to delocalized 
wavepackets at late times, while regions yielding substan- 
tial initial confinement, F > 0.9, correspond to localized 
wavepackets at late times. Additionally, we have seen 
good agreement between localized (delocalized) numeri- 
cal solutions and well (poorly) confined regions for other 
values of initial displacement that we have considered, 
\xc\ — 6 and 9.5. It should be noted that this model 
produces similar behavior to numerical and experi- 
mental [13] results reported previously for w > in the 
strongly interacting regime. While keeping u fixed, as the 
effect of anharmonicity is increased (in this case by mak- 
ing /? more positive, or in the previously reported cases 



by increasing the value of \xc\), the wavepacket becomes 
delocalized in the phase space, as seen in Fig. [T] 

We note that one of the interesting features of our 
model is the invariance of the appearance of localization 
under simultaneous sign changes of u and /3. This is re- 
flected in both the numerical results of the GP Equation, 
Fig. [TJd), and the topology of classical phase space. Fig. 
[2] The tendency to form a localized state is strongest 
when u and /3 have opposite signs: that is, when the 
trap oscillation frequency increases with action (/3 > 0) 
and interactions are attractive (u < 0) , or oscillation fre- 
quency decreases with action (/3 < 0) and interactions 
are repulsive (u > 0) . This latter case is an analog of the 
"negative mass instability" [28] that leads to bunching 
of a charged particle beam undergoing circular motion 
when the angular frequency decreases with energy. Par- 
ticles leading the bunch are repelled by the bunch, gain 
energy, lower their rotation rate, and fall back towards 
the bunch. By this process, a beam with a uniform dis- 
tribution of rotation phases will spontaneously form a 
bunch. Likewise, in a trap an initial wave function that 
is delocalized will spontaneously tend to bunch. 

Finally, we note that the interaction values considered 
here are accessible in experiments. One possible experi- 
mental realization is a lattice of one dimensional "tubes" 
formed by the interference of multiple lasers [29] . In our 
units, the required s-wave scattering length is related 
to the interaction parameter u as ag = uujQd/{2Nujj_), 
[18] , where N is the number of atoms per tube, uj± is the 
frequency in the transverse direction, and ground state 
width S = [h/{muJo)]^/^ ■ Taking uj±/{2Tr) - 49 KHz 
(which is reasonable with light of wavelength A = 1064 
nm and the atomic species ^^^K) ujo /{2Tr) — 2 KHz, 
N = 100, and u = 0.5 leads to ^ 0.036 nm. Re- 
cently, a group [nH [12] used the Feshbach resonance in 
•^^K to achieve as as small as 5.29 x 10"'^ nm . 
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